1 00:00:00,000 --> 00:00:02,770 NARRATOR: The following content provided by MIT 2 00:00:02,770 --> 00:00:06,400 OpenCourseWare under a Creative Commons license. 3 00:00:06,400 --> 00:00:08,505 Additional information about our license and 4 00:00:08,505 --> 00:00:11,800 MITOpenCourseWare in general is available at ocw.mit.edu. 5 00:00:11,800 --> 00:00:14,350 6 00:00:14,350 --> 00:00:21,180 PROFESSOR: So this is the second 086 lecture and it'll 7 00:00:21,180 --> 00:00:26,450 be my second and last lecture on solving ordinary 8 00:00:26,450 --> 00:00:33,270 differential equations -- u prime equal f of unt. 9 00:00:33,270 --> 00:00:41,030 So last time, we got Euler's method, both ordinary forward 10 00:00:41,030 --> 00:00:46,820 Euler, which is this method for the particular equation. 11 00:00:46,820 --> 00:00:51,440 So I'm using my model problem as always, u prime equal au. 12 00:00:51,440 --> 00:00:58,640 That's the model and I'm using small u for the true solution 13 00:00:58,640 --> 00:01:02,210 and capital u for the approximation. 14 00:01:02,210 --> 00:01:06,410 By whatever method we use, and right now, that method is 15 00:01:06,410 --> 00:01:09,050 forward Euler. 16 00:01:09,050 --> 00:01:14,570 We looked at the idea of stability, but now let me 17 00:01:14,570 --> 00:01:18,160 follow up right away with the point of stability. 18 00:01:18,160 --> 00:01:20,740 Where does stability come in? 19 00:01:20,740 --> 00:01:24,740 So I want to do that for this simplest of all methods, 20 00:01:24,740 --> 00:01:33,960 because you'll see how we get to an error equation and we'll 21 00:01:33,960 --> 00:01:35,780 get to a similar error equation 22 00:01:35,780 --> 00:01:38,710 for much better methods. 23 00:01:38,710 --> 00:01:44,320 So the accuracy of Euler is minimal, but the question of 24 00:01:44,320 --> 00:01:49,970 how stability enters and how we show that it converges will 25 00:01:49,970 --> 00:01:53,040 be the same also in partial differential equations, so 26 00:01:53,040 --> 00:01:55,670 this is the place to see. 27 00:01:55,670 --> 00:02:02,410 So I'll start with that, but then after that topic, will 28 00:02:02,410 --> 00:02:07,290 come better methods than Euler and the first step will be, of 29 00:02:07,290 --> 00:02:10,530 course, second order methods because Euler and backward 30 00:02:10,530 --> 00:02:14,290 Euler are only first order and we'll see -- actually, that'll 31 00:02:14,290 --> 00:02:15,300 be the point. 32 00:02:15,300 --> 00:02:19,470 Not only will show when they converge, why they converge, 33 00:02:19,470 --> 00:02:21,790 but also how fast they converge. 34 00:02:21,790 --> 00:02:23,460 What's the error? 35 00:02:23,460 --> 00:02:27,310 This estimate of the error is something that every code is 36 00:02:27,310 --> 00:02:31,490 going to somehow use internally to control the step 37 00:02:31,490 --> 00:02:33,530 size dela t. 38 00:02:33,530 --> 00:02:36,880 So I'm just going to go with the model problem, u prime 39 00:02:36,880 --> 00:02:40,260 equal au, the Euler method. 40 00:02:40,260 --> 00:02:45,830 So here is the evaluation -- au is f there, so it's really 41 00:02:45,830 --> 00:02:51,600 delta t multiplying f, which is just au in this easy case. 42 00:02:51,600 --> 00:02:56,700 Here's the true solution and then here is the new point. 43 00:02:56,700 --> 00:03:02,620 When I plug the true solution in to the difference equation, 44 00:03:02,620 --> 00:03:05,800 it, of course, doesn't satisfy exactly. 45 00:03:05,800 --> 00:03:09,310 There's some error term -- this, I would call the 46 00:03:09,310 --> 00:03:12,640 discreditzation error -- 47 00:03:12,640 --> 00:03:15,030 making the problem discrete. 48 00:03:15,030 --> 00:03:19,040 49 00:03:19,040 --> 00:03:24,840 So you could say the error at a single time step -- it's the 50 00:03:24,840 --> 00:03:28,510 new error that's introduced at time step m. 51 00:03:28,510 --> 00:03:32,370 52 00:03:32,370 --> 00:03:38,670 Error now is going to be -- e will be u, the true one, minus 53 00:03:38,670 --> 00:03:43,970 the approximate, so I just subtract to get an equation. 54 00:03:43,970 --> 00:03:49,960 For the error at time step n plus 1, that subtraction gives 55 00:03:49,960 --> 00:03:52,210 ena delta t. 56 00:03:52,210 --> 00:03:58,890 That subtraction gives en and this is the inhomogeneous 57 00:03:58,890 --> 00:04:02,620 term, you could say, the new error. 58 00:04:02,620 --> 00:04:05,320 So the question is -- 59 00:04:05,320 --> 00:04:10,340 first, the question is, what size is that, the local error? 60 00:04:10,340 --> 00:04:15,970 Then the key question is, when I combine and all these local 61 00:04:15,970 --> 00:04:20,460 errors, what's the global error? 62 00:04:20,460 --> 00:04:23,890 So first, the local error. 63 00:04:23,890 --> 00:04:28,310 You can see what this is, approximately, because the 64 00:04:28,310 --> 00:04:34,540 true solution would be e to the at and we all know that e 65 00:04:34,540 --> 00:04:39,760 to the at times u 0 -- 66 00:04:39,760 --> 00:04:40,990 e to the at -- 67 00:04:40,990 --> 00:04:44,930 there's a 1 plus an a -- 68 00:04:44,930 --> 00:04:47,470 from a single time step -- 69 00:04:47,470 --> 00:04:53,420 sorry, I should have said the true solution would grow by e 70 00:04:53,420 --> 00:04:56,700 to the a delta t over one step. 71 00:04:56,700 --> 00:05:00,950 It would multiply by e to the a delta t, but instead of 72 00:05:00,950 --> 00:05:04,030 multiplying by that exponential e to the a delta 73 00:05:04,030 --> 00:05:06,510 t, we're only multiplying by the first two 74 00:05:06,510 --> 00:05:08,940 terms in the series. 75 00:05:08,940 --> 00:05:17,960 So the error here is going to be something like a half delta 76 00:05:17,960 --> 00:05:23,020 t squared -- that's the key point where we see 77 00:05:23,020 --> 00:05:25,830 the local error -- 78 00:05:25,830 --> 00:05:29,220 times probably a squared -- 79 00:05:29,220 --> 00:05:33,300 80 00:05:33,300 --> 00:05:38,840 we needed an a squared and u. 81 00:05:38,840 --> 00:05:43,710 So a squared un, that would really be the second 82 00:05:43,710 --> 00:05:44,240 derivative. 83 00:05:44,240 --> 00:05:47,280 I'll put it in as second derivative, because that's 84 00:05:47,280 --> 00:05:53,390 what it would look like in a nonlinear case. 85 00:05:53,390 --> 00:05:57,980 So I'm working the linear case for simplicity, but the point 86 00:05:57,980 --> 00:06:03,360 is that this is the Euler local error. 87 00:06:03,360 --> 00:06:06,160 That's so typical, that you really want to look at that 88 00:06:06,160 --> 00:06:09,280 local error, see what it looks like. 89 00:06:09,280 --> 00:06:15,220 There is some constant and that has a certain effect on 90 00:06:15,220 --> 00:06:17,310 the quality of the method. 91 00:06:17,310 --> 00:06:20,240 There's a power of delta t and the question, what is that 92 00:06:20,240 --> 00:06:22,660 exponent has a big effect. 93 00:06:22,660 --> 00:06:29,430 So that exponent is the order plus 1. 94 00:06:29,430 --> 00:06:35,100 So for me, Euler's method has accuracy p equal 1. 95 00:06:35,100 --> 00:06:41,710 Over a single delta t step, the error is delta t squared, 96 00:06:41,710 --> 00:06:46,930 but we'll see that when I take 1 over delta t steps to get 97 00:06:46,930 --> 00:06:52,420 somewhere, to get to some fixed time like 1, then I have 98 00:06:52,420 --> 00:06:57,620 one over delta t of these things and the 2 comes back 99 00:06:57,620 --> 00:07:00,860 down to the 1, which I expect for Euler. 100 00:07:00,860 --> 00:07:07,550 And finally, this term is out of our control basically. 101 00:07:07,550 --> 00:07:11,680 That's coming from the differential equation. 102 00:07:11,680 --> 00:07:18,540 That would tell us how hard or easy the method is, the 103 00:07:18,540 --> 00:07:20,020 equation is. 104 00:07:20,020 --> 00:07:24,040 So if u double prime is real big, we're going to expect to 105 00:07:24,040 --> 00:07:28,470 have to reduce delta t to keep this error under control, so 106 00:07:28,470 --> 00:07:33,390 that's a typical term -- a constant delta t to a power of 107 00:07:33,390 --> 00:07:39,710 p plus 1 and the p plus first derivative, the second 108 00:07:39,710 --> 00:07:43,050 derivative matching that to the solution. 109 00:07:43,050 --> 00:07:44,810 So that's -- 110 00:07:44,810 --> 00:07:47,310 if I now put this -- 111 00:07:47,310 --> 00:07:49,990 112 00:07:49,990 --> 00:07:57,090 think of that as here as the sort of inhomogeneous term or 113 00:07:57,090 --> 00:08:00,680 the new source term. 114 00:08:00,680 --> 00:08:03,410 I just want to estimate the end now. 115 00:08:03,410 --> 00:08:06,950 So I've done the local part and now I'm interested in 116 00:08:06,950 --> 00:08:09,050 putting it all together. 117 00:08:09,050 --> 00:08:12,920 How do I look at e -- so I really want somehow the 118 00:08:12,920 --> 00:08:14,170 solution -- 119 00:08:14,170 --> 00:08:16,280 120 00:08:16,280 --> 00:08:21,240 e at some definite time n is what? 121 00:08:21,240 --> 00:08:28,030 I'm asking, really, for the solution to this system, and 122 00:08:28,030 --> 00:08:29,290 sort of in words -- 123 00:08:29,290 --> 00:08:33,380 124 00:08:33,380 --> 00:08:35,470 so that's a good -- 125 00:08:35,470 --> 00:08:39,600 a totally typical error equation. 126 00:08:39,600 --> 00:08:42,440 127 00:08:42,440 --> 00:08:47,800 I think the way to look at it in words is, at each step, 128 00:08:47,800 --> 00:08:53,440 some new term, new error is committed, and what happens to 129 00:08:53,440 --> 00:08:58,090 that error at later steps? 130 00:08:58,090 --> 00:09:02,660 At later steps, that error committed then gets multiplied 131 00:09:02,660 --> 00:09:06,990 by 1 plus a delta t at the next step and 1 plus a delta t 132 00:09:06,990 --> 00:09:08,780 at the following step. 133 00:09:08,780 --> 00:09:14,760 So I think we have something like after n steps, we would 134 00:09:14,760 --> 00:09:23,560 have 1 plus a delta t to the n times the error that we did -- 135 00:09:23,560 --> 00:09:24,810 our first step error. 136 00:09:24,810 --> 00:09:28,390 137 00:09:28,390 --> 00:09:31,010 So that's -- you see what's happening here? 138 00:09:31,010 --> 00:09:37,110 Whatever error we commit grows or decays according to this 139 00:09:37,110 --> 00:09:38,750 growth factor. 140 00:09:38,750 --> 00:09:41,980 Then of course, we made an error at step two and we made 141 00:09:41,980 --> 00:09:48,680 an error at step k, so at step k, it will have n minus k 142 00:09:48,680 --> 00:09:51,370 steps still to grow. 143 00:09:51,370 --> 00:10:04,010 This is de, the error de at step k and all the 144 00:10:04,010 --> 00:10:06,040 way through the -- 145 00:10:06,040 --> 00:10:10,440 I guess it would end with den. 146 00:10:10,440 --> 00:10:17,110 den would be the last thing that hadn't yet started to 147 00:10:17,110 --> 00:10:18,890 grow or decay. 148 00:10:18,890 --> 00:10:22,570 Ok so that's our formula. 149 00:10:22,570 --> 00:10:28,860 It looks a little messy, but it shows everything. 150 00:10:28,860 --> 00:10:31,920 It shows the crucial important of 151 00:10:31,920 --> 00:10:33,980 stability and how it enters. 152 00:10:33,980 --> 00:10:37,930 Stability controls -- so there is the error that we made way 153 00:10:37,930 --> 00:10:41,960 at the beginning and then this is the growth factor that 154 00:10:41,960 --> 00:10:49,870 happened n time, so you see that if 1 plus a delta t has 155 00:10:49,870 --> 00:10:54,890 magnitude bigger than 1, what's going to happen? 156 00:10:54,890 --> 00:11:01,960 That will explode and the computer will very quickly 157 00:11:01,960 --> 00:11:10,480 show totally out of bounds output. 158 00:11:10,480 --> 00:11:16,113 But if it's stable, if this thing is less than 1, less or 159 00:11:16,113 --> 00:11:21,160 equal to 1, then that stability means that this 160 00:11:21,160 --> 00:11:25,890 error is not growing and so we get -- 161 00:11:25,890 --> 00:11:29,840 so let me put, if we have this stability -- 162 00:11:29,840 --> 00:11:35,360 1 plus a delta t, smaller equal 1, then we get the bound 163 00:11:35,360 --> 00:11:35,980 that we want. 164 00:11:35,980 --> 00:11:38,500 So can you see what that looks like? 165 00:11:38,500 --> 00:11:41,500 166 00:11:41,500 --> 00:11:47,310 We've got n terms, right? 167 00:11:47,310 --> 00:11:50,510 To get to point n, we've made n -- 168 00:11:50,510 --> 00:11:52,600 capital n errrors -- 169 00:11:52,600 --> 00:11:57,140 and then each term, that gives a 1. 170 00:11:57,140 --> 00:12:04,980 That's the critical role of stability and the error is of 171 00:12:04,980 --> 00:12:10,270 this order, so I'm getting something like a half delta t 172 00:12:10,270 --> 00:12:12,640 squared over 2 times -- 173 00:12:12,640 --> 00:12:16,050 174 00:12:16,050 --> 00:12:18,350 let's say some maximum. 175 00:12:18,350 --> 00:12:24,390 I'll use u double prime for a maximum value, let's say, of 176 00:12:24,390 --> 00:12:27,580 the second derivative. 177 00:12:27,580 --> 00:12:34,000 That's a typical good satisfactory error estimate. 178 00:12:34,000 --> 00:12:36,060 So what's up here? 179 00:12:36,060 --> 00:12:41,140 So n times delta t times 1 delta t is the time that we 180 00:12:41,140 --> 00:12:41,840 reached, right? 181 00:12:41,840 --> 00:12:45,030 We took n steps, delta t each step. 182 00:12:45,030 --> 00:12:48,710 So this is the same as -- this is what I'm wanting and it's 183 00:12:48,710 --> 00:12:50,460 going to fit on this board -- 184 00:12:50,460 --> 00:12:54,070 there's a one half t -- 185 00:12:54,070 --> 00:13:01,000 I should make it t over 2, right? n delta t is some time 186 00:13:01,000 --> 00:13:03,460 -- capital t. 187 00:13:03,460 --> 00:13:04,830 We have the half. 188 00:13:04,830 --> 00:13:07,180 We still have 1 delta t. 189 00:13:07,180 --> 00:13:09,190 That's the key. 190 00:13:09,190 --> 00:13:14,650 We have this u double prime, which is fixed. 191 00:13:14,650 --> 00:13:18,410 Anyway, we have first order convergence. 192 00:13:18,410 --> 00:13:24,460 So the error after n steps goes down like delta t. 193 00:13:24,460 --> 00:13:27,540 194 00:13:27,540 --> 00:13:32,710 If I cut the time step in half, I cut the error. 195 00:13:32,710 --> 00:13:36,340 So that's not a great rate of convergence. 196 00:13:36,340 --> 00:13:38,530 That's, like, the minimal rate of convergence -- 197 00:13:38,530 --> 00:13:43,640 just 1 power delta t and that's why we call Euler a 198 00:13:43,640 --> 00:13:47,190 first order method. 199 00:13:47,190 --> 00:13:50,930 I think, I hope you've seen and can kind of -- 200 00:13:50,930 --> 00:13:56,760 because we'll -- this expressions for the solution 201 00:13:56,760 --> 00:14:00,110 is messy but meaningful. 202 00:14:00,110 --> 00:14:03,910 203 00:14:03,910 --> 00:14:10,850 That's a little bit of pure algebra or something that 204 00:14:10,850 --> 00:14:16,630 shows how stability is used, where it enters. 205 00:14:16,630 --> 00:14:22,070 So now I'm ready if you are to move to second order methods. 206 00:14:22,070 --> 00:14:24,690 207 00:14:24,690 --> 00:14:27,440 Now we're actually going to get methods that people would 208 00:14:27,440 --> 00:14:28,690 really use. 209 00:14:28,690 --> 00:14:32,340 210 00:14:32,340 --> 00:14:38,340 So all those methods are, I think, worth looking at. 211 00:14:38,340 --> 00:14:42,180 Can I say that the notes that are on the web -- you may have 212 00:14:42,180 --> 00:14:46,910 spotted chapter five, section one, for example, which is 213 00:14:46,910 --> 00:14:48,930 this material -- 214 00:14:48,930 --> 00:14:54,200 so I'm -- after the lecture, I always think back, what did I 215 00:14:54,200 --> 00:14:59,560 say and improve those notes, so possibly later this 216 00:14:59,560 --> 00:15:05,160 afternoon a revised, improved version of this material will 217 00:15:05,160 --> 00:15:10,060 go up and then next week starts the real thing, what I 218 00:15:10,060 --> 00:15:12,780 think of as the real thing, the wave equation and the heat 219 00:15:12,780 --> 00:15:16,020 equation -- partial differential equation -- but 220 00:15:16,020 --> 00:15:17,590 it's worth -- 221 00:15:17,590 --> 00:15:19,970 here we saw how stability works. 222 00:15:19,970 --> 00:15:28,140 Here we're going to see how accuracy can be increased and 223 00:15:28,140 --> 00:15:31,600 in different ways. 224 00:15:31,600 --> 00:15:34,120 In this way -- 225 00:15:34,120 --> 00:15:36,080 so what do I notice about the trapezoidal method? 226 00:15:36,080 --> 00:15:39,670 227 00:15:39,670 --> 00:15:43,140 Sometimes it's called Crank-Nickolson. 228 00:15:43,140 --> 00:15:47,010 In PDs, the same idea was proposed 229 00:15:47,010 --> 00:15:48,260 by Crank and Nickolson. 230 00:15:48,260 --> 00:15:50,840 231 00:15:50,840 --> 00:15:55,130 So what's their idea or what's the trapezoidal idea in the 232 00:15:55,130 --> 00:15:57,550 first place? 233 00:15:57,550 --> 00:16:00,240 Basically, we've centered -- 234 00:16:00,240 --> 00:16:04,710 we've recentered the method at what time position? 235 00:16:04,710 --> 00:16:07,280 236 00:16:07,280 --> 00:16:09,240 You can go ahead. 237 00:16:09,240 --> 00:16:10,770 So where -- 238 00:16:10,770 --> 00:16:16,080 instead of Euler, which was, like, started at time n or tn 239 00:16:16,080 --> 00:16:19,070 or backward Euler that was sort of focused on the new 240 00:16:19,070 --> 00:16:21,180 time tn plus 1? 241 00:16:21,180 --> 00:16:29,750 This combination is at time n plus a half delta t, right? 242 00:16:29,750 --> 00:16:36,090 Somehow we've symmetrized everything around the half and 243 00:16:36,090 --> 00:16:42,960 the point is that around the midpoint, the 244 00:16:42,960 --> 00:16:44,410 accuracy jumps up to 2. 245 00:16:44,410 --> 00:16:48,840 246 00:16:48,840 --> 00:16:51,350 Let me try always u prime equal au. 247 00:16:51,350 --> 00:16:54,990 248 00:16:54,990 --> 00:16:56,990 So f is au. 249 00:16:56,990 --> 00:17:01,030 Then what does this become? 250 00:17:01,030 --> 00:17:05,250 Can I, as you would have to do in coding it, separate the 251 00:17:05,250 --> 00:17:06,590 implicit part? 252 00:17:06,590 --> 00:17:07,930 This is implicit, of course. 253 00:17:07,930 --> 00:17:10,620 254 00:17:10,620 --> 00:17:17,650 That stands for this: fn plus 1 means f of the new unknown u 255 00:17:17,650 --> 00:17:21,060 at the new time. 256 00:17:21,060 --> 00:17:22,340 So it's implicit. 257 00:17:22,340 --> 00:17:26,760 258 00:17:26,760 --> 00:17:30,390 So let me bring that over to the left side, because it's 259 00:17:30,390 --> 00:17:31,510 sort of part of the unknown. 260 00:17:31,510 --> 00:17:37,490 I'll multiply through by delta t and I'll take f to be au. 261 00:17:37,490 --> 00:17:41,330 So I just want to rewrite that for our model problem. 262 00:17:41,330 --> 00:17:48,560 So it will be 1 minus a half when I bring that over minus a 263 00:17:48,560 --> 00:17:55,110 delta t over 2 un plus one. 264 00:17:55,110 --> 00:17:58,360 Is that what I get when I multiply through by delta t 265 00:17:58,360 --> 00:18:00,870 and bring the implicit part over and now I bring the 266 00:18:00,870 --> 00:18:06,060 explicit part to that side, so it's 1 plus a half 267 00:18:06,060 --> 00:18:08,870 delta t a delta tun. 268 00:18:08,870 --> 00:18:12,950 269 00:18:12,950 --> 00:18:17,000 So that's what we actually solve at every step. 270 00:18:17,000 --> 00:18:20,840 271 00:18:20,840 --> 00:18:25,210 So first of all, I see that it's pretty centered. 272 00:18:25,210 --> 00:18:27,530 I see that there's something has to be inverted. 273 00:18:27,530 --> 00:18:31,590 That's because the method's implicit. 274 00:18:31,590 --> 00:18:34,710 If this was a matrix, a system of equations, which is 275 00:18:34,710 --> 00:18:36,520 very typical -- 276 00:18:36,520 --> 00:18:40,230 could be a large system of equations with the matrix, 277 00:18:40,230 --> 00:18:47,200 capital a, then u would be a vector and we would be solving 278 00:18:47,200 --> 00:18:53,240 this system: matrix times vector equals known right hand 279 00:18:53,240 --> 00:18:55,650 side from time n. 280 00:18:55,650 --> 00:19:00,160 But here it's a scale or it's just a division, so what is 281 00:19:00,160 --> 00:19:03,700 stability going to depend on? 282 00:19:03,700 --> 00:19:05,680 What does stability require? 283 00:19:05,680 --> 00:19:08,960 At every step, we don't want to grow. 284 00:19:08,960 --> 00:19:15,230 So if we take a step -- that's 1 plus a over 2 delta t 285 00:19:15,230 --> 00:19:22,200 divided by 1 minus a over 2 delta t, that's the growth 286 00:19:22,200 --> 00:19:28,180 factor, the decay factor that multiplies un to give you n 287 00:19:28,180 --> 00:19:32,500 plus 1, so the stability question is, is this smaller 288 00:19:32,500 --> 00:19:34,600 or equal to 1? 289 00:19:34,600 --> 00:19:42,530 That's stability for this trapezoidal method. 290 00:19:42,530 --> 00:19:44,180 What's the story on that? 291 00:19:44,180 --> 00:19:48,580 Suppose a is very negative. 292 00:19:48,580 --> 00:19:51,880 That's what killed forward Euler, right? 293 00:19:51,880 --> 00:19:56,230 If a delta t -- that's where we look for instability. 294 00:19:56,230 --> 00:20:03,800 If a delta t becomes too negative, we lost stability in 295 00:20:03,800 --> 00:20:09,610 forward Euler and we'll lose stability in other, better 296 00:20:09,610 --> 00:20:10,890 methods too. 297 00:20:10,890 --> 00:20:17,140 But here, if a delta t is negative, that 298 00:20:17,140 --> 00:20:18,680 denominator is good. 299 00:20:18,680 --> 00:20:23,790 It's bigger than the numerator, right? 300 00:20:23,790 --> 00:20:25,520 So we have a -- 301 00:20:25,520 --> 00:20:26,630 you'll see it. 302 00:20:26,630 --> 00:20:31,940 If a delta t is less than 0, if it's very much less than 0, 303 00:20:31,940 --> 00:20:40,500 then this is approaching maybe minus 1, but it's below 1. 304 00:20:40,500 --> 00:20:44,860 So this message is a, stable. 305 00:20:44,860 --> 00:20:46,730 Absolutely stable. 306 00:20:46,730 --> 00:20:52,530 Stable for all a less than 0. 307 00:20:52,530 --> 00:20:55,140 We're fine. 308 00:20:55,140 --> 00:21:00,800 Delta t can be in principle as large as we like. 309 00:21:00,800 --> 00:21:05,850 So we might ask, why don't I just jump in one step to the 310 00:21:05,850 --> 00:21:09,630 time that I want to reach and not bother taking a lot of 311 00:21:09,630 --> 00:21:10,640 small steps? 312 00:21:10,640 --> 00:21:13,490 What's the -- 313 00:21:13,490 --> 00:21:15,270 in that little paradox? 314 00:21:15,270 --> 00:21:20,750 It would be stable, but what's no good. 315 00:21:20,750 --> 00:21:25,850 So if I jumped in one giant step, this thing would be 316 00:21:25,850 --> 00:21:27,630 about minus 1 or something. 317 00:21:27,630 --> 00:21:31,350 318 00:21:31,350 --> 00:21:34,660 So u here would be about minus 1 times u 0 and of course 319 00:21:34,660 --> 00:21:37,810 that's not the answer. 320 00:21:37,810 --> 00:21:44,030 So why do we still need a delta t? 321 00:21:44,030 --> 00:21:47,410 Because of the discretization error. 322 00:21:47,410 --> 00:21:57,780 Because this isn't the exact differential equation, so we 323 00:21:57,780 --> 00:21:59,940 still need -- 324 00:21:59,940 --> 00:22:03,600 we don't have problem with stability, but we still have 325 00:22:03,600 --> 00:22:13,080 to get the local errors down to whatever level we want. 326 00:22:13,080 --> 00:22:17,460 Now what are those local errors for this guy? 327 00:22:17,460 --> 00:22:25,040 What are the local errors for that message? 328 00:22:25,040 --> 00:22:25,530 Let's see. 329 00:22:25,530 --> 00:22:28,430 Can we actually see it? 330 00:22:28,430 --> 00:22:31,480 I mean, you know what I'm thinking, right? 331 00:22:31,480 --> 00:22:33,390 I'm thinking that the local error is a 332 00:22:33,390 --> 00:22:36,840 border delta t cubed. 333 00:22:36,840 --> 00:22:39,400 That's a second order method for me. delta t 334 00:22:39,400 --> 00:22:41,580 cubed at each step. 335 00:22:41,580 --> 00:22:45,860 1 over delta t steps, so delta t squared overall in that 336 00:22:45,860 --> 00:22:47,170 second order. 337 00:22:47,170 --> 00:22:49,780 But do we see why this is? 338 00:22:49,780 --> 00:22:52,920 I mean, instinctively, we say, it's centered. 339 00:22:52,920 --> 00:22:58,570 It should be, but let me expand this to see, are we -- 340 00:22:58,570 --> 00:23:02,450 341 00:23:02,450 --> 00:23:09,950 is e to the a delta t -- how close is it to 1 -- so it's 342 00:23:09,950 --> 00:23:17,690 approximately 1 plus a over 2 delta t divided by 1 minus a 343 00:23:17,690 --> 00:23:22,530 half a delta t. 344 00:23:22,530 --> 00:23:26,370 I just want to see that sure enough, this is second order. 345 00:23:26,370 --> 00:23:29,550 So of course I have a denominator here that I want 346 00:23:29,550 --> 00:23:30,730 to multiply. 347 00:23:30,730 --> 00:23:35,710 So that's 1 plus a delta t over 2, the explicit part and 348 00:23:35,710 --> 00:23:37,870 then everybody knows the 1 -- 349 00:23:37,870 --> 00:23:44,890 there are only two series to know, the truth is. 350 00:23:44,890 --> 00:23:48,170 We study infinite series a lot, but the two that 351 00:23:48,170 --> 00:23:52,630 everybody needs to know are the exponential series, which 352 00:23:52,630 --> 00:23:58,990 is for the left side and the geometric series for the right 353 00:23:58,990 --> 00:24:06,230 side, 1 over 1 minus x is 1 plus x plus x 354 00:24:06,230 --> 00:24:13,830 squared plus so on. 355 00:24:13,830 --> 00:24:16,980 356 00:24:16,980 --> 00:24:20,450 So that's the -- that this came from the denominator, 357 00:24:20,450 --> 00:24:24,180 getting it up in the numerator where I can multiply the other 358 00:24:24,180 --> 00:24:30,550 term and get 1 plus a delta t, which is good and what's the 359 00:24:30,550 --> 00:24:33,630 coefficient of a squared delta t squared? 360 00:24:33,630 --> 00:24:40,300 361 00:24:40,300 --> 00:24:41,630 Can you see what it is? 362 00:24:41,630 --> 00:24:45,530 If I do that multiplication, what does -- 363 00:24:45,530 --> 00:24:49,590 364 00:24:49,590 --> 00:24:50,330 you see it? 365 00:24:50,330 --> 00:24:54,980 Where do I get a delta t squared from this? 366 00:24:54,980 --> 00:25:00,980 I get one here, times 1/4 and I get 1 here times 1/4, so 367 00:25:00,980 --> 00:25:02,230 that's two 1/4s is 1/2. 368 00:25:02,230 --> 00:25:05,180 369 00:25:05,180 --> 00:25:07,040 So it's great, right? 370 00:25:07,040 --> 00:25:12,470 It got the next term correct in the exponential series. 371 00:25:12,470 --> 00:25:15,600 so that's why it's one order higher; because it got the 372 00:25:15,600 --> 00:25:18,720 right number there, but it won't get the right number at 373 00:25:18,720 --> 00:25:20,230 the next step. 374 00:25:20,230 --> 00:25:21,540 What's the -- 375 00:25:21,540 --> 00:25:25,000 I'd have to find out what the cube -- this would be a delta 376 00:25:25,000 --> 00:25:29,080 t over 2 cubed, so what will it -- 377 00:25:29,080 --> 00:25:31,250 the I'll get a wrong coefficient 378 00:25:31,250 --> 00:25:33,620 for a delta t cubed. 379 00:25:33,620 --> 00:25:36,310 I'll get -- 380 00:25:36,310 --> 00:25:40,950 that would give me 1/8 and this would give me another 381 00:25:40,950 --> 00:25:44,290 1/8, I think, so I think I'm getting two 1/8s out of that. 382 00:25:44,290 --> 00:25:47,290 I'm getting 1/4, which is wrong. 383 00:25:47,290 --> 00:25:49,660 What's right there? 384 00:25:49,660 --> 00:25:53,170 For the exponential series, it should be 385 00:25:53,170 --> 00:25:55,220 one over three factorial. 386 00:25:55,220 --> 00:25:56,820 It should be 1/6. 387 00:25:56,820 --> 00:26:03,460 So I got 1/4, not 1/6 and the error is the difference 388 00:26:03,460 --> 00:26:05,900 between them, which is 1/12. 389 00:26:05,900 --> 00:26:10,130 So 1/12 would be this number. 390 00:26:10,130 --> 00:26:17,470 The local error, so the local error, de, would be like the 391 00:26:17,470 --> 00:26:28,310 missing 1/12, the delta t to the third power now and the 392 00:26:28,310 --> 00:26:29,560 third derivative. 393 00:26:29,560 --> 00:26:32,340 394 00:26:32,340 --> 00:26:39,170 That's what would come out if I patiently went through the 395 00:26:39,170 --> 00:26:42,930 estimate, just to see, because we'll -- 396 00:26:42,930 --> 00:26:45,590 397 00:26:45,590 --> 00:26:49,990 in these weeks, we'll be creating difference methods 398 00:26:49,990 --> 00:26:53,700 for partial differential equations and we'll want to 399 00:26:53,700 --> 00:26:56,500 know, what's controlling the error? 400 00:26:56,500 --> 00:27:00,200 This is the kind of calculation, the beginning of 401 00:27:00,200 --> 00:27:03,370 a Taylor series that answers that question. 402 00:27:03,370 --> 00:27:09,890 So the main point is, local error delta t cubed, stable 403 00:27:09,890 --> 00:27:12,550 for every a delta t. 404 00:27:12,550 --> 00:27:19,170 So our method would move from Euler to trapezoidal and it 405 00:27:19,170 --> 00:27:22,770 would give us a final error e -- 406 00:27:22,770 --> 00:27:25,810 sorry, small e. 407 00:27:25,810 --> 00:27:27,680 I can say it without writing it. 408 00:27:27,680 --> 00:27:32,570 The error, small e, would be delta t to -- what power? 409 00:27:32,570 --> 00:27:34,680 So squared, delta t squared. 410 00:27:34,680 --> 00:27:36,690 Good. 411 00:27:36,690 --> 00:27:42,700 The same kind of reasoning -- now I want to ask about -- 412 00:27:42,700 --> 00:27:46,290 so that takes care of trapezoidal, which is a pretty 413 00:27:46,290 --> 00:27:47,980 good method. 414 00:27:47,980 --> 00:27:51,200 I think it's -- 415 00:27:51,200 --> 00:27:54,480 in METLAB, it would be ode 416 00:27:54,480 --> 00:27:58,220 something small t for trapezoid. 417 00:27:58,220 --> 00:27:59,670 Forgot -- maybe 1, 2. 418 00:27:59,670 --> 00:28:00,920 I should know. 419 00:28:00,920 --> 00:28:03,710 420 00:28:03,710 --> 00:28:04,960 It's used quite a bit. 421 00:28:04,960 --> 00:28:08,360 422 00:28:08,360 --> 00:28:12,250 Of course, we're going to get higher than second order, but 423 00:28:12,250 --> 00:28:16,850 second order is often in computational mathematics a 424 00:28:16,850 --> 00:28:18,310 good balance. 425 00:28:18,310 --> 00:28:22,350 Newton's method for solving nonlinear equation is second 426 00:28:22,350 --> 00:28:25,130 order and it doesn't pay to go to higher. 427 00:28:25,130 --> 00:28:26,910 You could invent a -- 428 00:28:26,910 --> 00:28:32,810 you could beat Newton by going, by inventing some third 429 00:28:32,810 --> 00:28:38,460 order method for solving a nonlinear equation, but in the 430 00:28:38,460 --> 00:28:41,240 end, Newton would be the favorite. 431 00:28:41,240 --> 00:28:43,600 So second order is pretty important. 432 00:28:43,600 --> 00:28:46,130 433 00:28:46,130 --> 00:28:50,390 So let me point to -- 434 00:28:50,390 --> 00:28:56,270 so that was implicit, because it involved fn plus 1, but 435 00:28:56,270 --> 00:29:04,250 Adams had another idea, which seems like a great idea. 436 00:29:04,250 --> 00:29:08,980 Instead of an implicit using fn plus 1, he used the 437 00:29:08,980 --> 00:29:14,590 previous value, fn minus 1, which we already know. 438 00:29:14,590 --> 00:29:20,240 We've already at that previous step substituted un 439 00:29:20,240 --> 00:29:23,520 minus 1 n to f. 440 00:29:23,520 --> 00:29:26,050 That might have been time consuming, but we sure had to 441 00:29:26,050 --> 00:29:31,480 do it once and we only have to do it once with Adams. 442 00:29:31,480 --> 00:29:33,690 We're using something we already know. 443 00:29:33,690 --> 00:29:37,540 This is the one new time that we need it, that we have to 444 00:29:37,540 --> 00:29:46,120 compute f and we get explicitly the new u. 445 00:29:46,120 --> 00:29:51,000 These numbers, 3/2 and minus 1/2 were chosen to make it 446 00:29:51,000 --> 00:29:56,510 second order and the notes check that, that with a series 447 00:29:56,510 --> 00:30:00,800 that it comes out to give us the correct second term. 448 00:30:00,800 --> 00:30:02,050 Good. 449 00:30:02,050 --> 00:30:04,450 450 00:30:04,450 --> 00:30:07,300 Let me mention backward differences. 451 00:30:07,300 --> 00:30:10,270 452 00:30:10,270 --> 00:30:17,730 Now we're just using 1 f and actually the implicit 1, so 453 00:30:17,730 --> 00:30:18,990 I'm making it implicit. 454 00:30:18,990 --> 00:30:23,390 So why am I interested in number 3 at all? 455 00:30:23,390 --> 00:30:26,030 Because it's going to have better stability. 456 00:30:26,030 --> 00:30:31,320 By being implicit and by choosing these numbers, 3/2, 457 00:30:31,320 --> 00:30:40,080 minus 1/2 and 1/2, I've upped the accuracy to 2, equal to 2 458 00:30:40,080 --> 00:30:43,470 and I've maintained high stability 459 00:30:43,470 --> 00:30:46,230 by making it implicit. 460 00:30:46,230 --> 00:30:53,700 So this one competes with this one for importance, for use. 461 00:30:53,700 --> 00:30:57,260 So you're really seeing three pretty good methods that can 462 00:30:57,260 --> 00:30:58,840 be written down. 463 00:30:58,840 --> 00:31:06,050 So these numbers were chosen to get that extra accuracy, 464 00:31:06,050 --> 00:31:09,590 which we didn't get from backward Euler. 465 00:31:09,590 --> 00:31:14,140 When those numbers were 1 and minus 1, backward Euler was 466 00:31:14,140 --> 00:31:15,080 only first order. 467 00:31:15,080 --> 00:31:22,490 Now this goes up to second order and you can guess that 468 00:31:22,490 --> 00:31:28,320 we can that every term we allow ourselves, we can choose 469 00:31:28,320 --> 00:31:34,880 our coefficients well, so that the order goes up by one more. 470 00:31:34,880 --> 00:31:45,150 Now if you look at that, you might say, why not keep -- 471 00:31:45,150 --> 00:31:47,870 you could come back to this. 472 00:31:47,870 --> 00:31:52,170 Why not combine the idea of backward differences, more 473 00:31:52,170 --> 00:31:58,150 left hand sides, with the Adams-Bashforth idea of more 474 00:31:58,150 --> 00:31:59,120 right hand sides? 475 00:31:59,120 --> 00:32:02,960 Suppose I come back to Adams-Bashforth and create an 476 00:32:02,960 --> 00:32:06,020 1806 method? 477 00:32:06,020 --> 00:32:09,390 That'll be third order, but still will 478 00:32:09,390 --> 00:32:11,090 only go back one step. 479 00:32:11,090 --> 00:32:15,320 So it will somehow use some combination like this that 480 00:32:15,320 --> 00:32:17,440 goes back -- 481 00:32:17,440 --> 00:32:21,385 I'm not writing this method down, for a good reason, of 482 00:32:21,385 --> 00:32:27,110 course, but it will use something like this on the 483 00:32:27,110 --> 00:32:34,360 left side and something like this on the right side and 484 00:32:34,360 --> 00:32:39,660 that gives us enough freedom, enough coefficients to choose 485 00:32:39,660 --> 00:32:44,890 -- two coefficients there, three there and enough that we 486 00:32:44,890 --> 00:32:48,620 can get third order accuracy. 487 00:32:48,620 --> 00:32:55,710 In other words, the natural idea is use both un minus 1 488 00:32:55,710 --> 00:33:00,320 and fn minus 1 in the formula to get that extra accuracy. 489 00:33:00,320 --> 00:33:04,740 So why isn't that the most famous method of all? 490 00:33:04,740 --> 00:33:08,430 Why isn't that even on the board? 491 00:33:08,430 --> 00:33:10,230 You can guess. 492 00:33:10,230 --> 00:33:16,810 It's violently unstable, so sadly, the most accurate 493 00:33:16,810 --> 00:33:22,980 methods, which use both an old value of u and an old value of 494 00:33:22,980 --> 00:33:27,630 f of u, use both of this and this -- 495 00:33:27,630 --> 00:33:32,030 the numbers that show up here are bad news. 496 00:33:32,030 --> 00:33:41,300 It's a fact of life and so that's why we go backwards. 497 00:33:41,300 --> 00:33:48,020 We have to make a choice: We go backwards with u or Adams 498 00:33:48,020 --> 00:33:56,350 goes backwards with f or we think of something way to hype 499 00:33:56,350 --> 00:34:02,020 up the accuracy even further within one step method. 500 00:34:02,020 --> 00:34:10,260 So the notes give a table, the beginning of a table. 501 00:34:10,260 --> 00:34:16,420 People could figure out the formulas for any order p, so 502 00:34:16,420 --> 00:34:19,615 the notes will give the beginning of a table for p 503 00:34:19,615 --> 00:34:23,410 equal 1, which is Euler, p equal 2, which is this. 504 00:34:23,410 --> 00:34:27,710 p equal 3 and p equal 4, which is Adams-Bashforth further 505 00:34:27,710 --> 00:34:36,350 back, backward differences further back and those tables 506 00:34:36,350 --> 00:34:39,200 show the constants. 507 00:34:39,200 --> 00:34:45,420 So they show the 1/12 or the 1/2 or whatever that is and 508 00:34:45,420 --> 00:34:49,840 they show the stability limit and of 509 00:34:49,840 --> 00:34:52,730 course, that's critical. 510 00:34:52,730 --> 00:34:55,870 The stability limit is much better for backward 511 00:34:55,870 --> 00:34:59,260 differences than it is for Adams-Bashforth. 512 00:34:59,260 --> 00:35:06,410 So actually, I only learned recently that Adams-Bashforth, 513 00:35:06,410 --> 00:35:11,580 which we all teach, which all books teach, I should say, 514 00:35:11,580 --> 00:35:15,400 isn't that much used way back -- maybe the astronomers might 515 00:35:15,400 --> 00:35:19,840 use it or they might use backward differences, which 516 00:35:19,840 --> 00:35:23,080 are more stable. 517 00:35:23,080 --> 00:35:26,240 So backward differences has an important role and then one 518 00:35:26,240 --> 00:35:28,680 step methods will have an important role. 519 00:35:28,680 --> 00:35:32,720 Backward differences are implicit, so those are great 520 00:35:32,720 --> 00:35:35,190 for stiff -- 521 00:35:35,190 --> 00:35:44,490 you turn that way for stiff equations and for nonstiff 522 00:35:44,490 --> 00:35:48,810 equations, let me show you what the workhorse method is 523 00:35:48,810 --> 00:35:51,000 in a moment. 524 00:35:51,000 --> 00:35:55,310 So I have a number 4 here, whose name I 525 00:35:55,310 --> 00:35:57,680 better put up here. 526 00:35:57,680 --> 00:36:04,880 Let me put up the name of method 4, which is two names: 527 00:36:04,880 --> 00:36:06,130 Runge-Kutta. 528 00:36:06,130 --> 00:36:08,260 529 00:36:08,260 --> 00:36:12,900 So that's a sort of one compound step. 530 00:36:12,900 --> 00:36:19,700 531 00:36:19,700 --> 00:36:21,840 So what do I mean by -- 532 00:36:21,840 --> 00:36:24,090 what's the point of -- 533 00:36:24,090 --> 00:36:28,470 this looked pretty efficient, but there are two aspects that 534 00:36:28,470 --> 00:36:33,100 we didn't really, I didn't really mention and on those 535 00:36:33,100 --> 00:36:36,650 two aspects, Runge-Kutta wins by being just a single 536 00:36:36,650 --> 00:36:39,330 self-contained step. 537 00:36:39,330 --> 00:36:46,190 So what are the drawbacks of a multistep method like 538 00:36:46,190 --> 00:36:47,950 this or like this? 539 00:36:47,950 --> 00:36:51,840 540 00:36:51,840 --> 00:36:57,130 You have to store the old value, no problem, but at the 541 00:36:57,130 --> 00:36:59,740 beginning, at t equals 0, what's the problem? 542 00:36:59,740 --> 00:37:02,420 543 00:37:02,420 --> 00:37:06,740 You don't know the old value, so you have to get started 544 00:37:06,740 --> 00:37:11,560 somehow with some separate -- you can't use this at t equals 545 00:37:11,560 --> 00:37:16,150 0, because it's calling for a value at minus delta t and you 546 00:37:16,150 --> 00:37:18,030 haven't got it. 547 00:37:18,030 --> 00:37:19,870 So you need -- 548 00:37:19,870 --> 00:37:21,580 it takes a separate method. 549 00:37:21,580 --> 00:37:23,950 You can deal with that, create -- use Runge-Kutta, for 550 00:37:23,950 --> 00:37:27,700 example, to get that first step. 551 00:37:27,700 --> 00:37:30,390 But then there's another problem with these backwards 552 00:37:30,390 --> 00:37:32,500 step methods, which again, can be resolved. 553 00:37:32,500 --> 00:37:36,760 So this is a live method. 554 00:37:36,760 --> 00:37:40,590 The other problem is, if you want to change delta t -- 555 00:37:40,590 --> 00:37:44,660 suppose you want to cut delta t in half, then this is 556 00:37:44,660 --> 00:37:49,450 looking for a value that's a half delta t back in time and 557 00:37:49,450 --> 00:37:54,220 you haven't got that either, but you've got values 1 delta 558 00:37:54,220 --> 00:38:00,680 t back and 2 delta t back and some interpolation process 559 00:38:00,680 --> 00:38:03,110 will produce that half value. 560 00:38:03,110 --> 00:38:04,260 So what am I saying? 561 00:38:04,260 --> 00:38:09,400 I'm just saying that getting started takes a special 562 00:38:09,400 --> 00:38:13,540 attention, changing delta t takes a little special 563 00:38:13,540 --> 00:38:15,710 attention, but still, it's all doable. 564 00:38:15,710 --> 00:38:21,340 So that this method, which I gave a name to last time -- 565 00:38:21,340 --> 00:38:23,270 I've forgotten what it was. 566 00:38:23,270 --> 00:38:36,960 I think it was something like ode 15 s for stiff is much 567 00:38:36,960 --> 00:38:40,150 used for stiff equations. 568 00:38:40,150 --> 00:38:42,410 So now I'm left -- 569 00:38:42,410 --> 00:38:47,040 let me not only write down Runge-Kutta's name, but the 570 00:38:47,040 --> 00:38:50,440 other -- so multistep I've spoken about. 571 00:38:50,440 --> 00:38:51,310 I just -- this is -- 572 00:38:51,310 --> 00:38:56,550 I'm just going to write this to remind myself to say one 573 00:38:56,550 --> 00:38:58,920 word about it -- 574 00:38:58,920 --> 00:39:00,170 dae. 575 00:39:00,170 --> 00:39:03,750 576 00:39:03,750 --> 00:39:05,460 Can I say one word even now? 577 00:39:05,460 --> 00:39:06,370 I'm sorry. 578 00:39:06,370 --> 00:39:09,030 Then I'll say one word later. 579 00:39:09,030 --> 00:39:14,570 So that's a differential algebraic equation. 580 00:39:14,570 --> 00:39:17,440 That's like a situation which comes up in chemical 581 00:39:17,440 --> 00:39:23,790 engineering and many other places, in which out of our n 582 00:39:23,790 --> 00:39:28,220 given equations, some of them are differential equations, 583 00:39:28,220 --> 00:39:33,370 but others are ordinary algebra equations, so there's 584 00:39:33,370 --> 00:39:38,270 no d by dt in them. 585 00:39:38,270 --> 00:39:40,590 Nevertheless, we have to solve and there is a d 586 00:39:40,590 --> 00:39:43,270 by dt in the others. 587 00:39:43,270 --> 00:39:47,830 So this is a whole little world of daes, which you may 588 00:39:47,830 --> 00:39:48,400 never meet. 589 00:39:48,400 --> 00:39:54,990 I have never personally met, but it has to be dealt with 590 00:39:54,990 --> 00:39:57,350 and there are codes -- 591 00:39:57,350 --> 00:40:03,900 I'll just mention the code dasil by [? pencil ?] 592 00:40:03,900 --> 00:40:05,830 I'll finish saying the one word. 593 00:40:05,830 --> 00:40:14,260 A model problem here would be some matrix times u prime 594 00:40:14,260 --> 00:40:23,460 equal f of u and t, so a sort of implicit differential 595 00:40:23,460 --> 00:40:26,100 equation, because I'm not giving u prime, I'm only 596 00:40:26,100 --> 00:40:28,700 giving mu prime. 597 00:40:28,700 --> 00:40:32,670 If m is singular -- 598 00:40:32,670 --> 00:40:37,970 otherwise I could multiply by m inverse and make it totally 599 00:40:37,970 --> 00:40:46,630 normal, but if m and sometimes t or m may depend on u and t, 600 00:40:46,630 --> 00:40:48,430 it could go singular. 601 00:40:48,430 --> 00:40:52,040 If it goes singular, then this system is degenerating and we 602 00:40:52,040 --> 00:40:53,900 have to think what to do. 603 00:40:53,900 --> 00:41:00,200 this dasil code would be one of the good ones that does. 604 00:41:00,200 --> 00:41:03,580 So that's sort of like my memory to myself to say 605 00:41:03,580 --> 00:41:09,940 something about the daes before completing this topic 606 00:41:09,940 --> 00:41:12,320 of ordinary differential equations. 607 00:41:12,320 --> 00:41:14,450 Now for Runge-Kutta. 608 00:41:14,450 --> 00:41:18,130 Let me take p equal to 2 first. 609 00:41:18,130 --> 00:41:22,130 The famous Runge-Kutta is the one at p 610 00:41:22,130 --> 00:41:23,430 with forth order accuracy. 611 00:41:23,430 --> 00:41:29,000 612 00:41:29,000 --> 00:41:35,990 That will take more of the little internal compounded 613 00:41:35,990 --> 00:41:39,570 steps than p equal to 2 takes, but p equal 614 00:41:39,570 --> 00:41:41,070 to 2 makes the point. 615 00:41:41,070 --> 00:41:45,300 616 00:41:45,300 --> 00:41:48,550 So this'll be simplified Runge-Kutta -- 617 00:41:48,550 --> 00:41:49,980 p equal to 2. 618 00:41:49,980 --> 00:41:55,030 It'll be un plus 1 minus un. 619 00:41:55,030 --> 00:42:01,830 You see, it's just one step, but what goes here? 620 00:42:01,830 --> 00:42:05,250 621 00:42:05,250 --> 00:42:06,260 Do I remember? 622 00:42:06,260 --> 00:42:07,930 The notes will tell me. 623 00:42:07,930 --> 00:42:13,010 I better -- since it's going to be kept forever, I better 624 00:42:13,010 --> 00:42:14,890 get it right. 625 00:42:14,890 --> 00:42:25,770 It involves fn and f of f, so Runge-Kutta, here we go. 626 00:42:25,770 --> 00:42:35,110 It has a half of fn, the usual Euler figure out this slope. 627 00:42:35,110 --> 00:42:43,690 But then if the other half, is the same f at -- 628 00:42:43,690 --> 00:42:48,210 I take Euler, so doing this part allowed me to taken a 629 00:42:48,210 --> 00:42:55,530 forward Euler and now I put that forward Euler step in as 630 00:42:55,530 --> 00:43:03,850 like a un plus 1 delta tfn at time -- 631 00:43:03,850 --> 00:43:04,830 what's the right time? 632 00:43:04,830 --> 00:43:07,040 Probably -- 633 00:43:07,040 --> 00:43:08,960 that's sort of -- 634 00:43:08,960 --> 00:43:14,240 this is like un plus 1, but we didn't have to do -- 635 00:43:14,240 --> 00:43:17,560 we're not implicit and so I presume that that's 636 00:43:17,560 --> 00:43:19,050 a time n plus 1. 637 00:43:19,050 --> 00:43:23,590 638 00:43:23,590 --> 00:43:28,210 So I've written in one line what the code would have to 639 00:43:28,210 --> 00:43:30,020 write in two lines. 640 00:43:30,020 --> 00:43:38,590 The first line of code would be compute f at the old un, 641 00:43:38,590 --> 00:43:44,400 then take an Euler step, multiply it by delta t and add 642 00:43:44,400 --> 00:43:51,450 to un, so this give something that's like un plus 1, but it 643 00:43:51,450 --> 00:43:55,420 didn't cost us anything. 644 00:43:55,420 --> 00:43:58,460 So can you compare that with trapezoidal? 645 00:43:58,460 --> 00:44:01,330 646 00:44:01,330 --> 00:44:05,780 You see the comparison with trapezoidal is trapezoidal has 647 00:44:05,780 --> 00:44:12,110 the 1/2 fn the same, but it has the 1/2 f n plus 1 at the 648 00:44:12,110 --> 00:44:18,940 new step involving the new un plus 1 where this 1 -- 649 00:44:18,940 --> 00:44:21,330 here is something like it, but something we 650 00:44:21,330 --> 00:44:23,600 already knew from Euler. 651 00:44:23,600 --> 00:44:27,120 So that's the natural idea. 652 00:44:27,120 --> 00:44:31,190 All these ideas actually would occur to all of us if we 653 00:44:31,190 --> 00:44:35,440 started thinking about these for a few weeks. 654 00:44:35,440 --> 00:44:40,010 B I guess I'm not planning to spend a few weeks on that, so 655 00:44:40,010 --> 00:44:46,610 I'm just jumping ahead to what people have thought of. 656 00:44:46,610 --> 00:44:51,680 So that's, you could say, two step Euler. 657 00:44:51,680 --> 00:44:56,200 There's two evaluations of f within a step. 658 00:44:56,200 --> 00:44:58,700 659 00:44:58,700 --> 00:45:00,130 So what's p equal 4? 660 00:45:00,130 --> 00:45:04,670 The famous Runge-Kutta is four of those. 661 00:45:04,670 --> 00:45:05,490 So you take -- 662 00:45:05,490 --> 00:45:06,500 there's an Euler step. 663 00:45:06,500 --> 00:45:08,300 That'll be step one. 664 00:45:08,300 --> 00:45:10,830 Then you'll stick that into something, 665 00:45:10,830 --> 00:45:13,890 that's step n plus 1/2. 666 00:45:13,890 --> 00:45:16,430 That'll be a second evaluation of f. 667 00:45:16,430 --> 00:45:18,450 That'll give you some output. 668 00:45:18,450 --> 00:45:23,560 You plug that in and that will actually, if it's correctly 669 00:45:23,560 --> 00:45:26,010 chosen, be another, better approximation of 670 00:45:26,010 --> 00:45:28,440 1/2, at n plus 1/2. 671 00:45:28,440 --> 00:45:30,810 You put that into the fourth one, which will give you 672 00:45:30,810 --> 00:45:34,180 something close to un plus 1 and then you 673 00:45:34,180 --> 00:45:35,900 take the right weights. 674 00:45:35,900 --> 00:45:37,150 Here they were 1/2 and 1/2. 675 00:45:37,150 --> 00:45:39,890 676 00:45:39,890 --> 00:45:42,330 Anyway, you can do it. 677 00:45:42,330 --> 00:45:46,360 The algebra begins to get messy beyond p equal 4, but 678 00:45:46,360 --> 00:45:48,730 you can go beyond p equal 4. 679 00:45:48,730 --> 00:45:53,400 I mean, there's books on Runge-Kutta methods because 680 00:45:53,400 --> 00:45:56,470 the algebra gets -- of these f of f of f. 681 00:45:56,470 --> 00:46:00,110 If there's anything that -- 682 00:46:00,110 --> 00:46:06,270 any construction in math that looks so easy, just take f of 683 00:46:06,270 --> 00:46:11,640 f of f of f of a function, of a number. 684 00:46:11,640 --> 00:46:14,510 That seems so simple, but actually it's extremely hard 685 00:46:14,510 --> 00:46:16,490 to understand what happens. 686 00:46:16,490 --> 00:46:22,270 Maybe you know that Mandelbrot's fractal sets and 687 00:46:22,270 --> 00:46:25,560 all these problems of chaos come exactly that way and 688 00:46:25,560 --> 00:46:29,940 they're extremely hard to see what's happening. 689 00:46:29,940 --> 00:46:34,900 When you do repeated composition, you would call 690 00:46:34,900 --> 00:46:37,760 that f of f of f. 691 00:46:37,760 --> 00:46:44,630 So Runge-Kutta stops at 4 and then there is actually an 692 00:46:44,630 --> 00:46:48,310 implicit Runge-Kutta that people work on, but I don't 693 00:46:48,310 --> 00:46:51,220 dare to begin to write that down and I won't even write 694 00:46:51,220 --> 00:46:55,760 the details of that one, which are in our notes and in every 695 00:46:55,760 --> 00:46:57,100 set of notes. 696 00:46:57,100 --> 00:46:58,370 But what's the point? 697 00:46:58,370 --> 00:47:01,920 698 00:47:01,920 --> 00:47:05,080 First, we get the accuracy, so what's the other thing you 699 00:47:05,080 --> 00:47:07,120 have to ask about? 700 00:47:07,120 --> 00:47:09,390 Stability. 701 00:47:09,390 --> 00:47:15,450 So what does our stability calculation give for 702 00:47:15,450 --> 00:47:16,700 Runge-Kutta? 703 00:47:16,700 --> 00:47:18,960 704 00:47:18,960 --> 00:47:21,500 The thing that you have to work with, actually, turns out 705 00:47:21,500 --> 00:47:23,150 to be quite nice. 706 00:47:23,150 --> 00:47:31,170 It's just the series for e to the at chopped off at the 707 00:47:31,170 --> 00:47:33,330 power, however far you're going. 708 00:47:33,330 --> 00:47:44,480 So this is Runge-Kutta stability, p equal to 2. 709 00:47:44,480 --> 00:47:51,250 The growth factor will be from -- if I just do this and I 710 00:47:51,250 --> 00:47:57,090 take the model problem, f equal au, you could -- 711 00:47:57,090 --> 00:48:01,170 if I took that out of there, the growth factor will be 1 712 00:48:01,170 --> 00:48:09,280 plus a delta t plus 1/2 a delta t squared. 713 00:48:09,280 --> 00:48:10,530 Stop. 714 00:48:10,530 --> 00:48:13,540 715 00:48:13,540 --> 00:48:17,560 That's what this would give for the model problem as the 716 00:48:17,560 --> 00:48:22,550 growth factor that multiplies each u to get the next u. 717 00:48:22,550 --> 00:48:26,890 What p equal 4 would give would be the 718 00:48:26,890 --> 00:48:31,990 same thing plus -- 719 00:48:31,990 --> 00:48:35,630 it's just the right -- it's just the exponential series 720 00:48:35,630 --> 00:48:43,370 chopped off at the fourth term. 721 00:48:43,370 --> 00:48:45,740 So of course, you see immediately, what's the 722 00:48:45,740 --> 00:48:47,840 discreditization error? 723 00:48:47,840 --> 00:48:49,420 It's still to t of the fifth. 724 00:48:49,420 --> 00:48:53,820 It's the first missing term that Runge-Kutta didn't 725 00:48:53,820 --> 00:49:00,330 capture, multiplied by 1 over 5 factorials, so you know 726 00:49:00,330 --> 00:49:03,740 right away that the constant is 1 over 120, 727 00:49:03,740 --> 00:49:07,280 which is pretty good. 728 00:49:07,280 --> 00:49:10,070 But the question is the stability. 729 00:49:10,070 --> 00:49:12,730 So what do we want to know here? 730 00:49:12,730 --> 00:49:15,350 We want to know -- 731 00:49:15,350 --> 00:49:18,170 so that means that a delta tb sum -- 732 00:49:18,170 --> 00:49:19,800 because that same number is coming up. 733 00:49:19,800 --> 00:49:21,810 Let me call it z, right? 734 00:49:21,810 --> 00:49:28,430 So the stability test for p equal to 2 is to look at the 735 00:49:28,430 --> 00:49:35,530 magnitude of 1 plus z plus 1/2 z squared and we want it to be 736 00:49:35,530 --> 00:49:43,370 less than or equal to 1. 737 00:49:43,370 --> 00:49:47,360 Now, the point is, this is something -- 738 00:49:47,360 --> 00:49:50,530 I hope you try this this weekend. 739 00:49:50,530 --> 00:49:54,440 Somehow get MATLAB to plot in the complex plane, 740 00:49:54,440 --> 00:49:56,180 so you really -- 741 00:49:56,180 --> 00:50:01,990 up there now we've just -- at some point z, some negative 742 00:50:01,990 --> 00:50:07,230 value of z -- it's going to hit 1 and then after that, 743 00:50:07,230 --> 00:50:11,040 it's going to be unstable. 744 00:50:11,040 --> 00:50:13,750 Of course, there's got to be a point where it's unstable. 745 00:50:13,750 --> 00:50:16,420 This is an explicit message. 746 00:50:16,420 --> 00:50:19,860 Runge-Kutta is explicit. 747 00:50:19,860 --> 00:50:23,200 So it can't -- if the problem is too stiff, we're going to 748 00:50:23,200 --> 00:50:26,840 be killed, but if it's an ordinary 749 00:50:26,840 --> 00:50:28,880 problem, this is a winner. 750 00:50:28,880 --> 00:50:35,600 So p equal 4, this is the code ode 45, which uses -- combines 751 00:50:35,600 --> 00:50:42,120 Runge-Kutta with 4 on a higher order formula to get an 752 00:50:42,120 --> 00:50:47,770 estimate for the error in that step. 753 00:50:47,770 --> 00:50:51,060 So I could make some comment on predictor, corrector 754 00:50:51,060 --> 00:50:54,100 methods, but I'll leave that in the notes. 755 00:50:54,100 --> 00:51:02,660 So what does that region look like? 756 00:51:02,660 --> 00:51:07,115 Then the other one is, what does the region 1 plus the z 757 00:51:07,115 --> 00:51:13,060 plus 1/2 c squared plus 1/6 c cubed plus 1/24 v to the 758 00:51:13,060 --> 00:51:18,760 fourth, less or equal 1 -- what's that region look like? 759 00:51:18,760 --> 00:51:24,680 I mean, if those regions are small, the method loses. 760 00:51:24,680 --> 00:51:28,580 If those regions are bigger than we get from 761 00:51:28,580 --> 00:51:31,990 Adams-Bashforth or backward differences, then the methods 762 00:51:31,990 --> 00:51:38,240 in there are good, successful. 763 00:51:38,240 --> 00:51:43,260 I've forgotten exactly what the picture is like there. 764 00:51:43,260 --> 00:51:45,000 Why am I in the complex plane? 765 00:51:45,000 --> 00:51:47,890 Because the number a -- 766 00:51:47,890 --> 00:51:50,430 z is a delta t. 767 00:51:50,430 --> 00:51:53,050 That number a can be complex. 768 00:51:53,050 --> 00:51:56,110 How can it be complex? 769 00:51:56,110 --> 00:51:59,520 I'm not really thinking that our model equation u prime 770 00:51:59,520 --> 00:52:05,040 equal au will be complex, but our model system would be u 771 00:52:05,040 --> 00:52:10,900 prime equals a matrix times u and then it's the Eigen values 772 00:52:10,900 --> 00:52:18,010 of that matrix that take the place of little a. 773 00:52:18,010 --> 00:52:22,370 So we have n different little as going for a system and the 774 00:52:22,370 --> 00:52:28,760 Eigen values of a totally real matrix can be nonreal, can be 775 00:52:28,760 --> 00:52:30,220 complex numbers. 776 00:52:30,220 --> 00:52:35,840 Let's just think when that happens and then I'll -- 777 00:52:35,840 --> 00:52:38,100 what does the figure looks like for p equals 4? 778 00:52:38,100 --> 00:52:41,910 I sort of remember that the figure even captures a little 779 00:52:41,910 --> 00:52:46,010 bit of that plane and then it goes -- 780 00:52:46,010 --> 00:52:48,430 why should I mess up the board with -- 781 00:52:48,430 --> 00:52:52,570 782 00:52:52,570 --> 00:52:54,440 it's a pretty good sized region. 783 00:52:54,440 --> 00:53:01,440 It's probably not shaped like that at all, but the point is, 784 00:53:01,440 --> 00:53:05,670 maybe this reaches out to about e, minus e. 785 00:53:05,670 --> 00:53:06,910 I think it does, somewhere around -- 786 00:53:06,910 --> 00:53:08,940 I don't know if that's an accurate -- 787 00:53:08,940 --> 00:53:11,960 right, reaches out around minus -- 788 00:53:11,960 --> 00:53:15,780 and reaches up to around 2. 789 00:53:15,780 --> 00:53:18,960 So in this last second, I was just going to say -- 790 00:53:18,960 --> 00:53:22,670 and we'll see it plenty of times -- where do we get 791 00:53:22,670 --> 00:53:27,810 complex -- where do we get nice, real, big systems, stiff 792 00:53:27,810 --> 00:53:36,360 systems with complex Eigen values in real calculations? 793 00:53:36,360 --> 00:53:40,560 By real by real calculations, I really mean pdes. 794 00:53:40,560 --> 00:53:44,850 So in a partial differential equation -- so this is just a 795 00:53:44,850 --> 00:53:49,350 throw away comment that we'll deal with again -- 796 00:53:49,350 --> 00:53:55,160 in partial differential equations, a uxx term or a uyy 797 00:53:55,160 --> 00:54:01,460 term that get multiply typically by diffusion, those 798 00:54:01,460 --> 00:54:07,920 are diffusion terms, those produce symmetric things, 799 00:54:07,920 --> 00:54:15,020 symmetric negative definitely, but convection, advection 800 00:54:15,020 --> 00:54:19,270 terms, velocity terms that get multiplied by some -- 801 00:54:19,270 --> 00:54:21,960 that are first derivatives -- 802 00:54:21,960 --> 00:54:25,890 those produce -- those are antisymmetric. 803 00:54:25,890 --> 00:54:28,940 Those produce these complex Eigen values that we 804 00:54:28,940 --> 00:54:30,630 have to deal with. 805 00:54:30,630 --> 00:54:36,480 So everybody wants to solve convection diffusion. 806 00:54:36,480 --> 00:54:40,170 Convection has these terms, diffusion has these. 807 00:54:40,170 --> 00:54:42,770 The coefficient may be very small. 808 00:54:42,770 --> 00:54:44,980 The hard -- that's the hard case, when the 809 00:54:44,980 --> 00:54:47,530 convection is serious. 810 00:54:47,530 --> 00:54:51,640 It's pushing us up the imaginary axis, where this is 811 00:54:51,640 --> 00:54:57,780 bringing us to the left and we're going to spend time 812 00:54:57,780 --> 00:55:00,350 thinking of good methods for that. 813 00:55:00,350 --> 00:55:02,570 So my homework -- 814 00:55:02,570 --> 00:55:05,280 would you like figure out this region and 815 00:55:05,280 --> 00:55:07,400 make a better pictures? 816 00:55:07,400 --> 00:55:14,370 Get MATLAB to do it and try any one of these methods on 817 00:55:14,370 --> 00:55:16,710 the model problem. 818 00:55:16,710 --> 00:55:18,470 Try a method or two. 819 00:55:18,470 --> 00:55:21,560 See whether the stability, the theoretical 820 00:55:21,560 --> 00:55:23,560 stability limit is serious. 821 00:55:23,560 --> 00:55:28,020 I mean, does the limit of a delta t -- 822 00:55:28,020 --> 00:55:29,270 can we -- 823 00:55:29,270 --> 00:55:31,450 824 00:55:31,450 --> 00:55:35,140 as I'm driving to MIT, I sometimes exceed the speed 825 00:55:35,140 --> 00:55:39,070 limit, the stability limit. 826 00:55:39,070 --> 00:55:43,420 Is that OK or not with finite differences? 827 00:55:43,420 --> 00:55:44,780 I'll see you Monday. 828 00:55:44,780 --> 00:55:45,800 Alright, thanks. 829 00:55:45,800 --> 00:55:47,480 Have a good weekend.